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Chapter 1 
Introduction 


1.1 Eigenvalue and Eigenvector Derivatives and their 

Applications 


Dynamic response and loads are an important consideration in the 
understanding and design of many physical systems. The analytical models 
fora wide range of these systems are governed by linear differential equations 
so that dynamic model analysis often consists of the solution of an eigenvalue 
problem. The eigenvalues and the eigenvectors of the system are fundamental 
quantities employed in determining the behavior of the system. Variations in 
system parameters lead to changes in the eigenvalues and the eigenvectors 
and hence in the response characteristics of the system. It is important to 


1 


know the magnitude of these variations, and this information is contained in 
the derivatives of the system eigenvalues and eigenvectors. Thus derivatives 
of eigenvalues and eigenvectors are of immense interest in several fields of 
physical sciences and engineering and much research effort has been 
expended in developing methods to calculate them. 

The applications of these derivatives (or synonymously, sensitivities) are 
manifold. Probably the most important applications are in the area of design 
optimization. System response sensitivities provide vital information in an 
optimization procedure and in general the cost of calculating derivatives is the 
dominant contributor to the total cost in an optimization procedure so that the 
efficient computation of eigenvalue and eigenvector derivatives is desirable. 
Derivatives can also be effectively used to approximate the eigenvalues and 
eigenvectors of a modified system and thus reduce the cost of reanalysis, 
substantially lowering the computational burden in optimization tasks. The 
derivatives are very useful even in non-automated design procedures because 
it is often not clear, from analysis alone, how to modify a design to improve 
or maintain the desirable properties. The derivatives identify design 
parameters that have the most or the least influence on the design process 
and thus ease the effort in design trend studies. 

Derivatives of eigenvalues and eigenvectors are particularly valuable in 
calculating the statistics of eigenvalue locations in stochastic dynamic 
systems. All physical systems are essentially subject to random environments 
and the effect of randomly changing environments is crucial for such systems 
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as missiles, spacecraft, airplanes, land vehicles, buildings and machinery. In 
addition, many system models do not have well-defined properties and it is 
frequently difficult to predict these properties (for example, stiffnesses) 
accurately[1-4]. The uncertainties in the system eigenvalues and eigenvectors 
are calculated from the estimated uncertainties in the properties of the system 
and the environment by using the derivatives of eigenvalues and eigenvectors. 

The application of derivatives is not restricted to design-oriented activities. 
Sensitivity analysis is also playing an increasing role in determining the 
analytical model itself In the areas of system identification and analytical 
model improvement using test results, sensitivity analysis is of growing 
importance. Much recent work in these fields is directly dependent on the 
calculation of eigensystem derivatives. 


1.2 Importance of higher order derivatives 


While in the past attention was mostly restricted to first order derivatives 
of eigenvalues, higher order derivatives are assuming a greater importance 
recently. It has been found in certain cases that second order derivatives are 
very effective in improving accuracy of approximations[5-13] and efficiency of 
design[7,1 1 ,12]. In almost all instances, eigenvalues are non-linear functions 
of system parameters and a second order approximation offers a wider range 
of applicability compared to the first order approximation. Intermediate 
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variables which may improve the quality of first order approximations are not 
generally available for eigenvalue approximations. Also, some optimization 
algorithms require second order derivatives, and first order derivatives of 
optimal solutions require second order derivatives of constraints[13]. The use 
of second derivatives can also greatly reduce the number of reanalyses 
required for the convergence of an optimization procedure[11,14]. Further, in 
certain optimization algorithms, second order approximations for eigenvalue 
constraints can drastically relax the move limits, thus achieving a nearly 
optimum trajectory, and can virtually eliminate the need for trial and error 
adjustment of move limits, thus improving the performance of the 
optimizer[14]. Looking at another aspect, in problems where instabilities are 
to be avoided, a first order calculation may completely fail to detect 
instabilities^]. References [15,16] also offer examples of the usefulness of 
second order derivatives. 


1.3 General Matrices 


The problem of calculating the derivatives of symmetric and hermitian 
eigenproblems is relatively simple and solution procedures are 
well-established, e.g. [17-21]. However, many physical problems give rise to 
non-self-adjoint formulations and thus lead to general matrices. An important 
example is aeroelastic stability which requires the solution of eigenproblems 
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with complex, general and fully populated matrices. General matrices are also 
obtained in damped structural systems and in network analysis and control 
system design where the eigenvalues are usually called poles. In the present 
study, the emphasis is on general matrices and the special properties of 
matrices, such as symmetry, are not considered. 

1.4 Approximate Updates to Eigenvalues 


Eigenvalue calculation for any but the smallest systems is an expensive 
process and is a major contributor to the computational expense of the typical 
dynamic analysis. In the design of structural systems, an iterative 
design/analysis process is performed until a satisfactory design is achieved. 
The cycle of the iterative process consists of an update of the structural 
design, a response analysis and calculation of updated responses and loads. 
When the process is not automated, a revision of the mathematical model may 
also be present. In typical dynamic analyses, each design cycle is an 
expensive process because the mathematical models are large and the time 
and effort required for analysis of an updated design are often prohibitive and 
greatly reduce the effectiveness of the design process. Apart from the 
computational effort, organizational effort can also be substantial. 

Appreciable cost savings can be realized if a quick evaluation of a change 
in system response resulting from the design changes is possible. An 
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efficient, even if approximate, evaluation of eigenvalues of a modified system 
is valuable in these applications. The value of an approximation depends on 
its efficiency as well as its accuracy in applications. A quick approximation 
that is valid in a very limited range of design space is of little use as it can 
severely reduce the global efficiency by requiring many more evaluations for 
convergence of the design process. 

Several researchers have worked on suitable approximations for 
eigenvalues of a modified design. However, in the past, attention seems to 
have been restricted to real symmetric systems which have eigenvalues in the 
real number field. It is one of the objectives of the present work to extend the 
techniques of approximation to general (complex non-hermitian) systems and 
perform a comparative analysis of the various techniques. 


1.5 Objectives of the Present Work 


The objectives of the present research are to: 

1. Review and perform a comparative analysis of the various methods 
available for calculating the derivatives of eigenvalues and eigenvectors 
of general matrices. 

2. Propose and evaluate some modifications to existing techniques. 
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understanding and design of many physical systems. The analytical models 
for a wide range of these systems are governed by linear differential equations 
so that dynamic model analysis often consists of the solution of an eigenvalue 
problem. The eigenvalues and the eigenvectors of the system are fundamental 
quantities employed in determining the behavior of the system. Variations in 
system parameters lead to changes in the eigenvalues and the eigenvectors 
and hence in the response characteristics of the system. It is important to 
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Know the magnitude of these variations, and this information is contained in 
the derivatives of the system eigenvalues and eigenvectors. Thus derivatives 
of eigenvalues and eigenvectors are of immense interest in several fields of 
physical sciences and engineering and much research effort has been 
expended in developing methods to calculate them. 

The applications of these derivatives (or synonymously, sensitivities) are 
manifold. Probably the most important applications are in the area of design 
optimization. System response sensitivities provide vital information in an 
optimization procedure and in general the cost of calculating derivatives is the 
dominant contributor to the total cost in an optimization procedure so that the 
efficient computation of eigenvalue and eigenvector derivatives is desirable. 
Derivatives can also be effectively used to approximate the eigenvalues and 
eigenvectors of a modified system and thus reduce the cost of reanalysis, 
substantially lowering the computational burden in optimization tasks. The 
derivatives are very useful even in non-automated design procedures because 
it is often not clear, from analysis alone, how to modify a design to improve 
or maintain the desirable properties. The derivatives identify design 
parameters that have the most or the least influence on the design process 
and thus ease the effort in design trend studies. 

Derivatives of eigenvalues and eigenvectors are particularly valuable in 
calculating the statistics of eigenvalue locations in stochastic dynamic 
systems. All physical systems are essentially subject to random environments 
and the effect of randomly changing environments is crucial for such systems 
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as missiles, spacecraft, airplanes, land vehicles, buildings and machinery. In 
addition, many system models do not have well-defined properties and it is 
frequently difficult to predict these properties (for example, stiffnesses) 
accurately[1-4]. The uncertainties in the system eigenvalues and eigenvectors 
are calculated from the estimated uncertainties in the properties of the system 
and the environment by using the derivatives of eigenvalues and eigenvectors. 

The application of derivatives is not restricted to design-oriented activities. 
Sensitivity analysis is also playing an increasing role in determining the 
analytical model itself In the areas of system identification and analytical 
model improvement using test results, sensitivity analysis is of growing 
importance. Much recent work in these fields is directly dependent on the 
calculation of eigensystem derivatives. 


1.2 Importance of higher order derivatives 


While in the past attention was mostly restricted to first order derivatives 
of eigenvalues, higher order derivatives are assuming a greater importance 
recently. It has been found in certain cases that second order derivatives are 
very effective in improving accuracy of approximations[5-13] and efficiency of 
design[7,11,12]. In almost all instances, eigenvalues are non-linear functions 
of system parameters and a second order approximation offers a wider range 
of applicability compared to the first order approximation. Intermediate 
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variables which may improve the quality of first order approximations are not 
generally available for eigenvalue approximations. Also, some optimization 
algorithms require second order derivatives, and first order derivatives of 
optimal solutions require second order derivatives of constraints[13]. The use 
of second derivatives can also greatly reduce the number of reanalyses 
required for the convergence of an optimization procedure[11,14]. Further, in 
certain optimization algorithms, second order approximations for eigenvalue 
constraints can drastically relax the move limits, thus achieving a nearly 
optimum trajectory, and can virtually eliminate the need for trial and error 
adjustment of move limits, thus improving the performance of the 
optimizer[14]. Looking at another aspect, in problems where instabilities are 
to be avoided, a first order calculation may completely fail to detect 
instabilities^]. References [15,16] also offer examples of the usefulness of 
second order derivatives. 


1.3 General Matrices 


The problem of calculating the derivatives of symmetric and hermitian 
eigenproblems is relatively simple and solution procedures are 
well-established, e.g. [17-21], However, many physical problems give rise to 
non-self-adjoint formulations and thus lead to general matrices. An important 
example is aeroelastic stability which requires the solution of eigenproblems 
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with complex, general and fully populated matrices. General matrices are also 
obtained in damped structural systems and in network analysis and control 
system design where the eigenvalues are usually called poles. In the present 
study, the emphasis is on general matrices and the special properties of 
matrices, such as symmetry, are not considered. 

1.4 Approximate Updates to Eigenvalues 


Eigenvalue calculation for any but the smallest systems is an expensive 
process and is a major contributor to the computational expense of the typical 
dynamic analysis. In the design of structural systems, an iterative 
design/analysis process is performed until a satisfactory design is achieved. 
The cycle of the iterative process consists of an update of the structural 
design, a response analysis and calculation of updated responses and loads. 
When the process is not automated, a revision of the mathematical model may 
also be present. In typical dynamic analyses, each design cycle is an 
expensive process because the mathematical models are large and the time 
and effort required for analysis of an updated design are often prohibitive and 
greatly reduce the effectiveness of the design process. Apart from the 
computational effort, organizational effort can also be substantial. 

Appreciable cost savings can be realized if a quick evaluation of a change 
in system response resulting from the design changes is possible. An 
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efficient, even if approximate, evaluation of eigenvalues of a modified system 
is valuable in these applications. The value of an approximation depends on 
its efficiency as well as its accuracy in applications. A quick approximation 
that is valid in a very limited range of design space is of little use as it can 
severely reduce the global efficiency by requiring many more evaluations for 
convergence of the design process. 

Several researchers have worked on suitable approximations for 
eigenvalues of a modified design. However, in the past, attention seems to 
have been restricted to real symmetric systems which have eigenvalues in the 
real number field. It is one of the objectives of the present work to extend the 
techniques of approximation to general (complex non-hermitian) systems and 
perform a comparative analysis of the various techniques. 


1.5 Objectives of the Present Work 


The objectives of the present research are to: 

1. Review and perform a comparative analysis of the various methods 
available for calculating the derivatives of eigenvalues and eigenvectors 
of general matrices. 

2. Propose and evaluate some modifications to existing techniques. 
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3. Formulate guidelines for selecting the most efficient computational 
algorithm for particular applications. 

4. Review the various approximations to eigenvalues for real symmetric 
systems and extend them to the case of complex general systems. 

5. Compare the various approximations in terms of efficiency and accuracy 
for some systems. 


1.6 Outline 


Chapter 2 reviews the various methods available for the calculation of 
derivatives of eigenvalues and eigenvectors for general matrices. The 
important consideration of normalization of eigenvectors of complex general 
matrices, which has not been adequately dealt with in the literature, is 
discussed. A new algorithm to calculate the derivatives of eigenvalues and 
eigenvectors simultaneously, based on a better normalizing condition, is 
described and important numerical aspects regarding the implementation of 
the algorithm are discussed, with consideration being given to sparse 
matrices. The various algorithms are classified as Adjoint or Direct. 

The efficiency considerations of the various algorithms are examined in 
Chapter 3. Operation counts are presented in terms of matrix size, number 
of design parameters, and the number of eigenvalues and eigenvectors of 
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interest. Actual CPU times are also presented for typical matrices for a range 
of parameters that influence the efficiency of the algorithms. 

Chapter 4 provides a survey of approximation methods proposed in the 
literature for real symmetric matrices, describing their special features. The 
approximation methods are extended to complex general matrices wherever 
feasible. Some approximation methods which do not seem to have been 
applied in the past are also presented. The approximation methods are 
classified on the basis of their theoretical origin. 

Numerical results from applying the proposed techniques of 
approximations are presented in Chapter 5. Operation counts are presented 
in terms of matrix size, number of design parameters, number of eigenvalues 
of interest and the number of times the approximation is to be performed. The 
approximation techniques are applied to typical matrices and random matrices 
and are evaluated in terms of their accuracy and efficiency. 

Chapter 6 contains the concluding remarks. General guidelines for 
selecting approximation methods and algorithms for calculation of eigenvalue 
and eigenvector derivatives are summarized. This chapter also contains 
remarks about the limitations of the present work and recommendations for 
further research. 
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Chapter 2 

Derivatives of Eigenvalues and Eigenvectors for 

General Matrices 


2.1 Problem Definition 


The matrix eigenproblem is defined as follows: 

Au W = \ {k) u {k) (2.1.1) 

and the corresponding adjoint problem is 

v W7 A = l {k V k)T ( 2 . 1 . 2 ) 
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where A is a general complex matrix of order n and A.W, uW and vM are the 
k -th eigenvalue and right and left eigenvectors respectively. A superscript T 
denotes the transpose. 

(The adjoint problem is defined by some authors in an alternative form as 

,<*>A - nV 

where superscript * denotes a conjugate-transpose. However, the notation of 
eq. (2.1.2) is more popular in the literature on structural dynamics). 

The eigenvalues and eigenvectors are complex and do not necessarily 
occur in complex-conjugate pairs. All eigenvalues are assumed to be distinct. 

The matrix A and hence, A.W, uM and vW are functions of design 
parameter vector p with individual parameters denoted by Greek subscripts, 
e.g. p a . Derivatives with respect to p H are denoted by the subscript , a e.g., 
= A _ . All the design variables are assumed to be real. 

dp a 

The well-known biorthogonality properties of the eigenvectors are given 
by 

v«V> = 0 iff / # j (2.1.3) 


and 


v (,)7_ Au^ = 0 iff / j 


(2.1.4) 
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Note that, the left hand side of eq. (2.1.3) is not an inner product as usually 
understood, since vW and/or may be complex vectors. Note also that the 
left eigenvectors of A are the right eigenvectors of A 7 and vice versa. 


2.2 Normalization of eigenvectors 


The eigenvectors uM and vW are not completely defined by eqs. (2.1.1) 
and (2.1.2). A normalization condition has to be imposed to obtain unique 
eigenvectors. For brevity, let us consider only the normalization of the right 
eigenvector. A normalizing condition frequently imposed in the self-adjoint 
case is the following: 


u wy*> = i 


( 2 . 2 . 1 ) 


However, it is not always possible to use eq. (2.2.1) for non-self-adjoint 
problems as can equal zero or a very small number causing numerical 


difficulties. 


his is true even if the matrix A is real, as shown by the example 
0 -1 


matrix A = 


Unfortunately, considerable confusion exists in the 


H oj 

literature regarding this point and several authors arbitrarily adopted 
eq. (2.2.1) as a normalizing condition for non-self-adjoint problems, 
e.g.[1 1,12, 15, 22-25], In this respect, the formulations of these references are 
not rigorous for general matrices. 

One possible way to avoid the above difficulty is to replace eq.(2.2.1) by 
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u W‘u<*> = 1 


( 2 . 2 . 2 ) 


where superscript * denotes a conjugate-transpose. Eq. (2.2.2) is not prone to 
the difficulties of eq. (2.2.1) because uW’uM is always guaranteed to be 
non-zero. But, eq. (2.2.2) is not a complete normalizing condition as it does not 
render the eigenvector unique. If u satisfies eq. (2.2.2), then w = ue' c , where 
/ = V ~ 1 and c is an arbitrary real number, also satisfies eq. (2.2.2). Despite 
this limitation, eq. (2.2.2) can be used satisfactorily in certain 
formulations[26,27]. 

Another normalization condition, inspired by the biorthogonality properly 
of the left and right eigenvectors, is 


v wy*) = ! 


(2.2.3) 


Eq.(2.2.3) also does not render the eigenvectors unique, since a pair 
uW, vM can be replaced by cuM,(1/c)vW, where c is an arbitrary non-zero 
complex number, and still satisfy eq. (2.2.3). Again, this is not necessarily a 
severe restriction for calculation of the derivatives of eigenvectors^, 16, 28-30], 
It must, however, be emphasized that if the eigenvector is not unique, nor is 
its derivative. 

The normalization condition 

uff = 1 (2.2.4) 
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is very attractive because it renders the eigenvectors unique and at the same 
time, the index m can be chosen easily to avoid ill-conditioning. Apparently, 
only Nelson[31] used this normalizing condition in obtaining the derivatives 
of eigenvectors. The index, m, may be chosen such that 

luff I = mpx | tiff | (2.2.5) 

Another choice for m, used by Nelson[31], is 
luff | Ivffl = max luff | |vff| (2.2.6) 


The nature of uncertainty of the derivative of the eigenvector is of some 
interest. Without a normalizing condition, an eigenvector is uncertain to the 
extent of a non-zero constant multiplier. The derivative of an eigenvector is 
uncertain to the extent of an additive multiple of that,eigenvector. To show 
this, let uM be an eigenvector so that wM = cuM is also an eigenvector. 
Then, if p a is a design parameter, 


<W*> <W‘>) a„<*> . „,(*) 

— = r = c— + dll' 

°Pn d Pa d Pa 


(2.2.7) 


dc 

The quantity d = — — is not zero since the quantity c is not really a 

dp a 

constant, but is a function of the nature of the normalization criterion. In 
practice, the constant d depends on the way the eigenvectors uM and wW are 
normalized. 
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2.3 Methods of Calculation 


The various methods of calculating the derivatives of eigenvalues and 
eigenvectors can be divided into three categories: 

1. Adjoint Methods, which use both the right and the left eigenvectors. 

2. Direct Methods, which use only the right eigenvectors. 

3. Iterative Methods, which use an iterative algorithm that converges to the 
required derivatives. 


2.3.1 Adjoint Methods 


The first expressions for the derivatives of eigenvalues of a general matrix 
seem to have been derived by Lancaster[32). Considering only a single 
parameter, Lancaster obtained the following expressions for the first and 
second derivatives of an eigenvalue: 


X 


(k) = 

, (I 


v<*> r A w u<*> 


v wy*> 


(2.3.1) 
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(2.3.2) 


*(*) _ 
A .aa _ 


Jk ) T a U W 

v A . aa u n 

+ 21 

vwy/o j= i 

j*k 


(yW r A a u«>) (v^A a u<*>) 

(x<*) - ( V w v>) 


Eq. (2.3.1) can be obtained in the following manner. Differentiate eq. (2.1.1) 
with respect to the parameter p a to obtain 


A a u (,[) + Au« - X< k a V*> + 


(2.3.3) 


Premultiplying both sides by vW T , we get 


v wr A a u<*> + v<*>Wa “ + v<«W« 


(2.3.4) 


The last terms in the expressions on both sides are equal due to eq. 
(2.1.2), so that 


r<*> r A a u<*> - X.W V W 7 uW 


(2.3.5) 


Eq. (2.3.1) follows immediately. Eq. (2.3.2) is derived in its more general form 
later. An expression corresponding to eq. (2.3.1) for a non-linear eigenvalue 
problem 

A(A,)u W = 0 (2.3.6) 

was obtained by Pedersen and Seyranian[33] in a similar manner as 
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(2.3.7) 


V W 7 A u W 

{ k ) • « 


Ka = 


v<*> t A x u 


(*) 


To obtain the second derivatives of eigenvalues, the first derivatives of left 
and right eigenvectors are calculated either explicitly[5,12,16,26,30] as in eq. 
(2.3.16) or implicitly[12,15,32] as in eq. (2.3.18). Since the eigenvalues are 
assumed to be distinct, the set of eigenvectors forms a basis for the n-space 
and the first derivatives of eigenvectors can be expressed in terms of the 
eigenvectors as 


u 


(*) _ 


n 

1 c 

j= 1 


kja 


u 


(/) 


and 


.W = 


£ * -« 

j = i 




(2.3.8) 


Now, the calculation of the first derivatives of eigenvectors reduces to the 
evaluation of the coefficients c k j a and d k j a . 

Premultiplying eq. (2.3.3) by v^ T , where j ^ k , we get 




(2.3.9) 


Now, substituting the expansions of eq. (2.3.8) and using eqs. (2.1.1) and (2.1.2) 


and the bi-orthogonality property of eq. (2.1.3), we obtain 


c kja ~ 


v 0 t A u w 

", (I U 




k *j 


(2.3.10) 
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Proceeding in a similar manner after differentiating eq. (2.1.2) with respect to 
the parameter p a , we obtain 


^kjn 


v<« r A a uW 

(A.W - Jt w )v« r U W 


k*j 


(2.3.11) 


The above expressions for the coefficients c^ ra and t/^ (l were obtained by 
Rogers[29j. 

It can be observed that 


^kja 


yWV*) 

Cy *° v (/) V 


(2.3.12) 


Reddy[34] derived an equivalent expression for the response derivative by 
casting the derivative as the solution of a forced response problem for the 
same system. 

Note that, in view of eq. (2.2.7), the coefficients c kka and d kkn in eq. (2.3.8) 
are arbitrary and depend on the normalization of the eigenvectors. For 
example, if eq. (2.2.4) is used to normalize the right eigenvectors, then 


c kka 


V r /|VJ 
. ^ . { 'kja u m 

7=1 

j*k 


(2.3.13) 


and if eq. (2.2.3) is used to normalized the left eigenvectors, then 


^ kka c 


kka 


(2.3.14) 
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It has been proposed[31,35,36] that the eigenvector derivative be 
approximated by using less than the full set of eigenvectors in the expansion 
of eq. (2.3.8) so that the evaluation of eigenvector derivative by Adjoint method 
could become cheaper. This variant of Adjoint method has received mixed 
reports in the literature[31,35]. The quality of such an approximation is difficult 
to assess beforehand and the selection of the number of eigenvectors to be 
retained in the expansion is problem dependent. It is not considered in this 
work because a meaningful comparison with other methods cannot be easily 
be made. However, this consideration should not be ignored while 
implementing the sensitivity calculations for particular problems. 

The expressions for the second derivatives of eigenvalues were obtained 
by Plaut and Huseyin[30]. For the sake of simplicity in expressions, let us 
assume, without loss of generality, that the right and left eigenvectors are 
normalized as in eq. (2.2.3). Eq. (2.3.1) can then be written as 

X!*> - v< * )r A a u ( *> (2.3.15) 


Differentiating with respect to a parameter uncorrelated to the parameter 
p a , we obtain 


iW =Jk)T . (*) 

, o(t v , «|t u 


+ 




+ 


{k)T. (k) 

v ; 0 A , a u 


(2.3.16) 


which can be equivalently written, without involving the derivative of the left 
eigenvector, as 
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- v'‘ lr A, „ (i u W + » ( ‘ ,r (A „ - X + «<*>„ - X f lW» (2.3.17) 


Eq. (2.3.16) can be rewritten using eqs. (2.3.10) and (2.3.11) as 


*!«ft = + y 5 1 (^’ - *•'">«**<% + CjypC/jy,.) 

j*k 


.(*) 


i(*) - i</K 


(2.3.18) 


Crossley and Porter[5,28] derived similar expressions for derivatives with 
respect to the elements of the matrix. Expression for the N-ih order diagonal 
derivative was derived by Elrazaz and Sinha[9] and it is 


+ + 


N y (k) 


d n x 


<K 


N 


= v Wr 


°Pa 


^r 1 ^ 

/V(/V - 1) /~ 2 a d 2 u {k) 

. N-2 -2 

°Pa oPa 

-'A *W-1 (*) 

<;A o u' ' 


2! 


+ 


... + N 


0p„ X'.N-'l 

ra op a 

N(N — 1) d N ~\ (k) f) 2 u W 
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... - /V 


^r 2 

<ix w gw-yw 

% 5 P «-’ 


(2.3.19) 


Morgan[37] developed a different computational approach for the 
derivative of an eigenvalue without requiring the eigenvectors explicitly. His 
expression is equivalent to 
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(2.3.20) 


, W = 

\ a 


trace of ^[adj(A - X^I)]a a J- 
trace of adj(A - X^l) 


The corresponding expression for derivatives with respect to matrix elements 
was derived by Nicholson[38], 

It can however be shown that[39] 

adj(A - X {k) \) = t k u {k V k)T (2.3.21) 


where t k is a constant and that[40] 

trace of -£[adj(A - X^I)]a = t k v ^ T A (1 i/^ 
trace of adj(A - X^i) = t k v^ T 


(2.3.22) 


Thus, in the computation of adj(A - XM|), both right and left eigenvectors 
are implicitly computed, in view of eq. (2.3.21). Eqs. (2.3.22) also show that 
Morgan's eq. (2.3.20) is equivalent to Lancaster's eq. (2.3.1). Woodcock[41] 
also obtained formulas involving the adjoint matrix for the first and second 
derivatives of eigenvalues. An operation count shows that calculation of the 
adjoint matrix is several times more expensive than the explicit calculation of 
right and left eigenvectors so that Lancaster's formula is preferable to 
formulas requiring the adjoint matrix. This conclusion is also supported by 
sample computations[42]. In addition, although eq. (2.3.20) was used 
satisfactorily for small problems[43,44], numerical difficulties were reported for 
reasonably large problems[45]. Woodcock's formula for the second derivative 
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of the eigenvalue requires a partial derivative of the adjoint matrix and this is 
so complicated that Woodcock himself recommends the finite difference 
method. Formulas due to Morgan and Woodcock are not therefore considered 
in the following. 

In calculating the derivatives by Adjoint Methods, i.e., using eqs. (2.3.1), 
(2.3.8)-(2.3.18), 

• the first derivative of an eigenvalue requires the corresponding right and 
left eigenvectors. 

• the first derivative of an eigenvector requires all the left and right 
eigenvectors. 

• the second derivative of an eigenvalue requires the corresponding right 
and left eigenvectors and their first derivatives. 

2.3.2 Direct Methods 

The second category comprises methods that evaluate the derivatives 
using only the right eigenproblem. Direct methods typically involve either the 
evaluation of the characteristic polynomial or the solution of a system of linear 
simultaneous equations without requiring all the left and right eigenvectors. 
Methods requiring the evaluation of the characteristic polynomial and the 
derivative of the determinant[45,46] are 0(n 5 ) processes while other methods 
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considered here are at most 0(n 3 ) processes. In addition, the determination 
of the characteristic polynomial is, in general, an unsatisfactory process with 
respect to numerical stability, even when all the eigenvalues are 
well-conditioned[47]. While numerically stable algorithms have been 
proposed for evaluation of the characteristic polynomial[48], the computational 
expense still seems to be formidable. Hence, we do not consider these 
methods. Methods requiring the solution of a system of equations have the 
particularly attractive feature that the coefficient matrix needs to be factored 
only once for each eigenvalue regardless of the number of parameters and the 
order of the derivatives required. Thus, they are very useful in applications 
where higher order derivatives are required. 

The earliest method in this class is due to Garg[22] who obtained the first 
derivatives of the eigenvalue and the eigenvector by solving two systems of 
(n + 1) equations each in the real domain, without requiring any left 
eigenvectors. However, his formulation involves several matrix 
multiplications. Rudisill[23] proposed a scheme in which only the 

corresponding left and right eigenvectors are required to calculate the first 
derivative of the eigenvalue and the eigenvector. This was refined by Rudisiil 
and Chu[24] to avoid calculating the left eigenvectors altogether. Solution of 
a system of only (n + 1) equations is required (though in the complex domain) 
to obtain the first derivatives of eigenvalue as well as eigenvector. Extension 
to higher order derivatives is straightforward. Cardani and Mantegazza[25] 
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proposed solution methods of the same formulation for sparse matrices and 
extended it to the quadratic eigenproblem. 

One weakness that is common to all the above formulations that do not 
require left eigenvectors[22-25] is that they rely on the normalization condition 
given by eq. (2.2.1), which is unreliable as discussed earlier. 

Nelson[31] circumvented this difficulty by using the normalizing conditions 

v(* ,T u w = 1 and u$ = 1 (2.3.23) 


However, the formulation of Rudisill and Chu is superior to .No!; 


formulation in that it does not require any left eigenvectors. 

In this work, we propose a variation of the Rudisill and Chu formulation 
which does not rely on the questionable normalizing condition of eq. (2.2.1) 
and at the same time requires no left eigenvectors. 

Differentiating eq. (2.1.1), we get 


Au ( *J + A 


U W _ X W U W + 

a u , a 


xV 


(2.3.24) 


which can be rewritten in partitioned matrix form as 


[a - x<*>i : 



A, a u 


(*) 


(2.3.25) 


Now, we impose the normalizing condition of eq. (2.2.4). Differentiation of 
eq. (2.2.4) yields, 
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(2.3.26) 


Because of eq. (2.3.26), the m-th column of the coefficient matrix in eq. 
(2.3.25) can be deleted. Eq. (2.3.26) also reduces the number of unknowns by 
one so that eq. (2.3.25) is now a system of n equations in n unknowns. Eq. 
(2.3.25) is rewritten as 

By 1 = r (2.3.27) 


where 


B = [A - X {k) \ 


_ (kU 

u J/rMh column deleted 


Vi = 



with /77-th clement deleted 


r= “ A ,a u 


(k) 


(2.3.28) 


To get second derivatives, differentiate (2.3.24) with respect to and get, 


/a _ iW hll W _ _ _ 

(A A l)U a p U A_ aB - 


ap 


A .aP u 


(k) _ /a _ >Wn„(k) 

l A ,a p 


-(Ap - 


(2.3.29) 


or, in partitioned matrix form, 
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(2.3.30) 


u {k) 


— A u W — (A — l Wn.-W 

, a(i u ' , a ' U ,H 


", a|l 


(A. |l - 


Following the same reasoning as before, eq. (2.3.30) is written as 


By 2 = s 


(2.3.31) 


where 


V2 = 



with m-th element deleted 


s - - A. all" 1 * 1 - (A, a - X»»I)U^ - (A | t - A<>u[*> (2.3.32) 

Note that, if is a distinct eigenvalue of A and if ujj {> * 0, then the matrix 
A is of rank (n — 1) and the m-th column that is deleted is linearly dependent 
on the other columns. Hence the matrix B is non-singular. The matrix B will 
also be well-conditioned if uffl is the largest component in the eigenvector 
uM and the matrix A is itself not ill-conditioned. The vectors y., and y 2 can be 
obtained by standard solution methods. If the matrix A is banded or if the 
derivatives of both right and left eigenvectors are required, it may be more 
efficient to use a partitioning scheme as described below. 
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2.3.3 Modification of Direct Method for Banded Matrices 


Equations (2.3.27) and (2.3.31) can be written as 


(A A/ ^!)m-th column deleted 1 *, a m-th row deleted ^ r (2.3.33) 


Let uW be normalized so that uffl = = constant 

Eq. (2.3.33) is a system of n equations. Writing the m-th equation 
separately, we have, if the superscript ( k ) is omitted for notational 
convenience, 


Cx , « “\« x =t 


a m x , a \ a u mO 


(2.3.35) 


where 


C — (A - Xl) m _. 


m-th row and column deleted 


x , a u ,a m-th row deleted 


x u m-th row deleted 


r m-th row deleted 
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aj„ = m - th row of A with the m-th column deleted 
From (2.3.35), 

\ a u mO = a m x , a ~ r m 
From (2.3.34), 

x a = C _1 (\ a x + t) 

Eliminating x a , we have 

— r m 

i — m m 

, (Z T 

u m0 — x 

where 

— [C ] a m 

Proceeding in a similar manner for the left eigenvector, 

y,„“[cYV.„y + *,) 

where 

y,a = v , a m-th row deleted 
y ~ v m-th row deleted 


(2.3.36) 


(2.3.37) 


(2.3.38) 


(2.3.39) 
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t/ ( r /)m-th row deleted 


r t being the appropriate right hand side. 

Thus the following procedure can be used to obtain the derivatives A a and 

u .a- 

1. Form a LU decomposition of the matrix C. 

2. Solve b m = [C T ] _ 1 a m by forward substitution. 

3. Calculate A. a from (2.3.38). 

4. Calculate x n from (2.3.37) by backward substitution. 

5. Expand x „to u n setting u m a = 0. 

If the derivatives v a of the left eigenvectors are also required, only three 
further steps are needed. 

6. Calculate y a from (2.3.39) by forward substitution. 

7. Expand y to v (1 setting v m a = 0. 

8. Normalize v n appropriately depending on the normalization of v. For 
example, to obtain the derivative of the left eigenvector that satisfies the 
normalization condition of eq. (2.2.3), subtract (v r u , a + v , T a u)v . 
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The matrix C needs to be factored only once. Also, the matrix C retains 
the bandedness characteristics of the original matrix A, so that advantage can 
be taken of it. Furthermore, higher order derivatives can be obtained by 
merely substituting an appropriate right hand side vector, r. However, higher 
order derivatives can suffer in accuracy because of accumulated round-off 
error. 

The conditioning of matrix C needs some comment. Note that C is 
obtained from the singular matrix (A - A,W|) by deleting both the row and 
column corresponding to index m. Hence, for matrix C to be non-singular, one 
must make sure that the m-th row is linearly dependent on the other rows as 
well as that the m - th column is linearly dependent on the other columns. In 
other words, C is non-singular iff u$ * 0 and # 0 . If vffl is very small 
compared to the largest element in vW, steps 2 and 4 in the above procedure 
will give inaccurate results even if u$> is the largest element in u^. In 
general, it is not possible to make a good choice for m without the knowledge 
of the left eigenvector. Since the calculation of left eigenvector using forward 
substitution in an inverse iteration scheme is cheap (as explained later in 
Section 3.1), it is suggested that the left eigenvector be calculated and the 
index m be chosen as in eq. (2.2.6). This is the same criterion used by 
Nelson[31] and will assure as well-conditioned a matrix C as possible. 

In summary, we note that, in calculating derivatives by Direct Method, 

• left eigenvectors are not used. 


29 


• a complete solution of the eigenvalue problem is not required, if the 
derivatives of only a few of the eigenvalues and eigenvectors are sought. 
This is in contrast to Adjoint Method which requires all the left and right 
eigenvectors to calculate the first derivative of any eigenvector. 

• calculation of any derivative requires the solution of a system of linear 
equations. 

• only one matrix factorization needs to be performed for all orders of 
derivatives of an eigenvalue and its corresponding right and left 
eigenvectors. 

2.3.4 Iterative Methods 

Andrew[27] proposed an iterative algorithm to calculate the first 
derivatives of eigenvalues and eigenvectors. This algorithm is a refined and 
generalized version of the iterative scheme developed by Rudisill and Chu[24], 
Except for the dominant eigenvalue, the convergence of this algorithm seems 
to be very much dependent on the choice of the initial values for the 
derivatives. To be efficient for non-hermitian matrices, this iterative method 
requires a complex eigenvalue shifting strategy which is not easy to 
implement. Hence this method is not considered. 
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Chapter 3 


Efficiency Considerations in Calculating Derivatives 


In order to establish criteria for the selection of the most efficient algorithm 
for calculating the derivatives of eigenvalues and eigenvectors in a given 
application, we compare the operation counts and actual CPU times required 
by Adjoint and Direct methods. The decision as to which algorithm is best is 
necessarily problem-dependent. The comparison is, however, described in 
terms of three variables that can usually be ascribed to a given problem and 
which significantly influence the decision. These variables are 

1. the size of the matrix n 

2. the number of design parameters m 

3. the number of eigenvalues of interest I. 
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3.1 Operation Counts 


To start with, let us consider the operation counts (multiplications and 
divisions only) for the adjoint methods given by eqs. (2.3.1),(2.3.8)-(2.3.18) and 
the direct methods given by eqs.(2.3.27)-(2.3.32). They are summarized in 
Table 1. It should be noted that the operation counts represent an estimate 
of the actual number of operations performed by a solution routine and include 
only the most significant terms. The actual number of operations will vary 
slightly depending on programming details. The effect of the sparsity of the 
matrix derivative A a is modeled in the operation counts by the parameter k, 
defined such that the number of operations in evaluating the product A a u is 
equal to K/7 2 (that is, k = 1 corresponds to a full A a ). 

The eigenvalues are calculated using the EISPACK subroutine package[49] 
by first reducing the matrix to upper hessenberg form using unitary similarity 
transformations and then applying the QR algorithm. The number of 
operations and the CPU time for calculating the eigenvalues is not relevant in 
evaluating the methods to calculate the derivatives. The operation count for 
eigenvalue computation is given only for comparison. 

The right eigenvectors are calculated by inverse iteration on the same 
upper hessenberg matrix used for calculating the eigenvalues and are back 
transformed using standard subroutines in the package EISPACK. The 
corresponding operation count is given in Table 1. The inverse iteration 
algorithm is an extremely powerful method for computing eigenvectors and is 


32 



much superior in accuracy as well as speed of convergence to the common 
alternative algorithms based on the solution of homogeneous equations or 
direct iteration. Algorithms based on the solution of homogeneous equations 
are limited in their accuracy by the accuracy of the eigenvalue and those 
based on direct iteration are limited in their convergence, particularly for 
eigenvectors not corresponding to either the largest or the smallest 
eigenvalue. 

For the calculation of left eigenvectors, it is important to note that there is 
no need to repeat the process with the transposed matrix. The left 
eigenvectors are obtained cheaply using forward substitution in place of 
backward substitution in the inverse iteration process. There is also no need 
to repeat the matrix factorization. A subroutine was written to calculate the 
right and left eigenvectors in this manner and the corresponding operation 
count is given in Table 1. 

Table 1 gives the operation count of evaluating the individual steps. To 
obtain the number of operations involved in evaluating the derivatives, we 
must add the operation counts for all the steps required in the calculations. 
These counts are given in the following discussion. 
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Table 1. Operation Counts 


Eigenvalues and Eigenvectors 


Operation Operation Count 


Evaluation of eigenvalues 
Evaluation of right eigenvectors 

Evaluation of left eigenvectors 


8 n 3 to 10n 3 
H2n 2 ) 

4 n2 ) 


Adjoint Methods 


Operation 


Operation Count 


Evaluation of eq. (2.3.1) 
Evaluation of eq. (2.3.8), 
(2.3.10), (2.3.11) 

Evaluation of eq. (2.3.18) 


Imn 2 k 
lmn 2 { ic + 2) 



n 2 K 


Direct Methods 


Operation Operation Count 


LU decomposition of matrix B 
Formulation and solution of eq. (2.3.27) 

Formulation and solution of eq. (2.3.31) 


/m$ 2 ( k + 1) 

I n 2 ( 3k + 1) 
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3.2 CPU Time Statistics 


In the following tables, computational cost for the calculation of the first 
and second derivatives of eigensystems are compared for matrices of order 
20, 40 and 60. The CPU time statistics are obtained on the IBM 3084 computer 
using the VS-FORTRAN compiler with no compiler optimization. The 
correlation between operation counts and CPU times is shown in Tables 2 and 
3. The ratio of operation count(OC) and CPU time for various operations, 
tabulated in Tables 2 and 3, is about 10 5 operations per CPU second with a 
variablity of 27 percent. 

The typical matrices are generated for the dynamic stability aeroelastic 
analysis of a compressor stage rotor with mistuned blades. The geometric and 
structural parameters of the rotor and formulation and method of analysis are 
the same as those of NASA Test Rotor 12 described in reference[50] except 
that the number of blades and the torsional frequencies are varied. The 
torsional frequency values are selected randomly from a population of mean 
1.0 and standard deviation 0.01. The standard deviations of the actual 
samples are slightly different. 
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Table 2. Correlation between Operation Count(OC) and CPU time 


Calculation of right eigenvectors 


n 

1 


OC/CPU seconds (xIO 4 ) 

60 

60 


8.6 

60 

10 


8.3 


Calculation of right and left eigenvectors 

n 

1 


OC/CPU seconds (xIO 4 ) 

60 

60 


8.5 

60 

10 


9.2 



Evaluation of eq. 

(2.3.1) 

n 

1 

m 

OC/CPU seconds (xIO 4 ) 

60 

60 

10 

10.7 

60 

60 

5 

10.7 

60 

10 

10 

10.7 

60 

10 

5 

10.7 


Evaluation of eq. (2.3.8), (2.3.10), (2.3.11) 

n 

1 

m 

OC/CPU seconds (xIO 4 ) 

60 

60 

10 

12.0 

60 

60 

5 

11.9 

60 

10 

10 

11.9 

60 

10 

5 

11.9 
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Table 3. Correlation between Operation Count(OC) and CPU time(Contd.) 


Evaluation of eq. (2.3.18) 

n I m OC/CPU seconds (xIO 4 ) 

60 60 10 9.1 

60 10 5 7.4 

Decomposition of matrix B 

n I OC/CPU seconds (xIO 4 ) 

60 60 8.5 

60 10 9.2 

Evaluation of eq. (2.3.27) 

n I m OC/CPU seconds (xIO 4 ) 

60 60 10 12.9 

60 10 5 13.0 




3.3 Calculation of First derivatives of Eigenvalues only 


Operation Count 


Adjoint Method 
Direct Method 


l(-^-n 2 + k mn 2 ) 
l[~ + (k + \)mn 2 ] 


It is clear from the operation count that the Adjoint Method, which is an 
0(n 2 ) process, is superior to Direct Method, an 0(/7 3 ) process, for large n. The 
number of design variables and the number of eigenvalues of interest have no 
bearing on this conclusion. As the order of the matrix increases, the direct 
method becomes more expensive. For example, for 5 design variables and 
10 eigenvalues of interest, the CPU time for the Direct method is 2.3 times 
more expensive than for the Adjoint Method for n = 20, and for n = 60, the 
ratio is 3.0. 
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3.4 Calculation of First derivatives of Eigenvalues and 
Eigenvectors 


Operation Count 

Adjoint Method 
Direct Method 


-^-n 3 + lmn 2 ( k + 2) 

. n 

-y- + lmn 2 { k + 1) 


When the derivatives of both eigenvalues and right eigenvectors are 
required, the choice of method is dependent on the values of / and m. When 
very few eigenvalues are of interest, the Direct method is cheaper. When 
many eigenvalues are of interest, the Direct method is more expensive than 
the Adjoint method. However, this effect of the number of eigenvalues of 
interest is less significant when the number of design variables is large. As 
the number of design variables increases, the direct method becomes more 
competitive, even when all eigenvalues are of interest. For a 60 x 60 full 
(k = 1) matrix, this is illustrated in Figure 1 on page 40. 

The operation count shows that the computation by adjoint method of 
eigenvector derivative, which is necessary for the second derivative of 
eigenvalue, is an 0(/? 3 ) process and is more expensive than the computation 
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of the eigenvector itself which is an 0(n 2 ) process using the procedure 
described in Section 3.1. This fact is significant as some authors have stated 
the opposite[6,7]. 


41 



3.5 • Calculation of First and Second derivatives of 
Eigenvalues only 


Operation Count 

Adjoint Method 
Direct Method 
Direct-Adjoint Method 


-^-n 3 + (k + 1 )mn 3 4- / j n 2 ic 
/ (2 )« 2(3k + 1) 

/-— -I- / j n 2 K + lmn 2 (2K + 1) 


The Direct-Adjoint Method denotes the calculation of the eigenvector 
derivatives by the Direct method and the eigenvalue second derivatives by the 
Adjoint Method. The third term in the operation count for the Direct-Adjoint 
Method is significant only when m is small. From the operation count, it is 
seen that the Direct-Adjoint Method is always cheaper than the Direct Method. 
Hence, the choice lies between the Adjoint Method and the Direct-Adjoint 
method. Here, considerations similar to those of the last section hold and the 
choice of method depends on the values of / and m. When few eigenvalues 
are of interest, the Direct-Adjoint method is cheaper. When many eigenvalues 
are of interest, the Adjoint method is superior. But this advantage of Adjoint 
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CPU I IMP. (SiruNDS) 



Figure 2. CPU Times for calculation of second derivatives of eigenvalues for a 60 x 60 matrix 
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Chapter 4 

Approximate Eigenvalues of Modified Systems 


4.1 Introduction 


The eigenvalue problem to be solved is 


Au<*> = *<*>«.<*> 


( 4 . 1 . 1 ) 


and 


v<‘> r A - X<*V‘» T 


( 4 . 1 . 2 ) 


where A is a general complex matrix of order n and uM and vW are the 
k -th eigenvalue and right and left eigenvectors respectively. 
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The matrix A differs by a small amount from a nominal matrix A 0 . The 
eigenvalues X 0 and the right and left eigenvectors u 0 and v 0 are taken to be 
known. 

A = A 0 + AA (4.1.3) 


A „W - iW«W 
A 0 U 0 - ^0 u 0 


(4.1.4) 


and 




i(*U*)r 

A o v o 


(4.1.5) 


Throughout this chapter, we will also assume that the left eigenvectors are 
normalized such that 


Jk)T (k) _ 1 
v 0 u 0 - i 


(4.1.6) 


The eigenvalues X and the eigenvectors u and v can be written as 

= A#> + AX W 

u W = ujf } + Au W (4.1.7) 

v {k) = v<*> + Av (A) 


Let p ra0 be the vector of nominal design variables and Ap ra be a vector of 
perturbations from the nominal design variables so that 
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A 0 “ A (P«o) 


(4.1.8) 


and 


A = A(p a0 + Ap a ) (4.1.9) 

Approximate quantities are denoted by the subscript a. For example, an 
approximation for the eigenvalue XW is denoted by X^ h K All derivatives are 
evaluated at the nominal design. 

An exact relation exists between AX, AA, Au and Av as follows. 

> W + vW T AAu<*> + 4*>(v<f> 7 Au W + Av WT u 0 ) 

+ vjf )T AAAu (fc) + Av W7 A 0 Au W 

lk) + Av (fc)7 AAuj, fc) + Av (/<)7 AAAu {k) 

X {k) (4.1.10) 

1 + Av<*> r u 0 + VgAu + AvW r AuW 

The object is to obtain without solving a full eigenvalue problem. X.W 
can be obtained exactly using eq. (4.1.10), if Au and Av are Known. Since the 
exact values Au and Av cannot be obtained without solving a full eigenvalue 
problem, various approximations can be formed based on the above 
expression. 

The approximations are broadly classified into 

1. Derivative based approximations 

2. Rayleigh-Quotient based approximations 
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3. Trace-theorem based approximations 


** * 

4. Others 


4.2 Derivative Based Approximations 


Derivative based approximations are of special importance in optimization 
problems because first derivatives are required anyway in most optimization 
algorithms. 

The most common of the derivative based approximations are based 
directly on truncated Taylor series. We will consider Linear and Quadratic 
approximations in this category. The enormous cost associated with the 
computation of any higher derivatives with multiple design variables 
effectively precludes the possibility of using higher order approximations 
based on the Taylor series. 

4.2.1 Linear Approximation(LIN) 

This is the simplest approximation to be considered in this work. Linear 
Approximation is obtained by truncating the Taylor series expansion for the 
eigenvalue after two terms. 
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(4.2.1) 


(*) _ >(*) , 
'■a — A 0 ^ 


m 

I 


dX 


(k) 


a = 1 dp a 


-&Pa 


A linear approximation is usually inadequate in terms of accuracy because 
eigenvalues are often highly non-linear functions of design variables. The 
linear approximation will be referred to herein as LIN approximation. 


4.2.2 Quadratic Approximation(QUAD) 


The Quadratic approximation is obtained by truncating the Taylor series 
expansion for the eigenvalue after three terms. 


y{k) _ *(*) 
A a — An 


+ 


(IX 


(k) 


a=1 (7p a 


~ A Pa + T 


m 

I 


d 2 X {k) 


2 d= i|{ = i <lp a cpfi 


Ap (I Ap (i 


(4.2.2) 


The quadratic approximation can be quite expensive for large orders of the 
matrix or large number of design variables. Miura and Schmit[14] used a 
simplified form of the quadratic approximation and found that, considering 
global efficiency, the higher cost of the quadratic approximation can 
sometimes offset the higher accuracy in an optimization problem. The 
quadratic approximation will be referred to herein as QUAD approximation. 

Because of these efficiency considerations, attempts were made to 
improve the accuracy of the linear approximations through the use of 
intermediate variables with respect to which the eigenvalues may be nearly 
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linear. However, Miura and Schmit[14] concluded that such intermediate 
variables cannot be found for a general structural problem. 

The accuracy of linear approximations can also be improved in another 
fashion. This is by introducing non-linearities without, however, introducing 
the second derivatives, which are expensive to calculate. 


4.2.3 Conservative Approximation 


In optimization applications, it is often desired to have a conservative 
approximation. For eigenvalue problems, this usually means underestimating 
the eigenvalues. Starnes and Haftka[51] proposed a hybrid approximation, 
which is a combination of a linear approximation and a reciprocal 
approximation (linear in 1/p (I ) such that it is the most conservative 
combination of the two. Since this approximation is only applicable to real 
quantities, it is applied to the real part and/or the imaginary part of the 
eigenvalue as required in particular applications. 


l(« - iW + 
A. ^ — /vn » 


m 

i 

a=i 


e. 


dx 


W 


d Pa 


A P« 


(4.2.3) 


where 
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Sk™ 

1 if — r ^ 0 

d Pa 

PaO dX {k) 

if — > 0 


Pa 


°P a 


Even though this approximation is not, in general, more accurate than the 
linear approximation, it is popular because it is more conservative than the 
linear approximation and it is also convex. 


4.2.4 Generalized Inverse Power Approximation 


Non-linearities can also be introduced into the linear approximation in a 
more direct manner, as in the Generalized Inverse Power approximation, 
described by Prasad[52,53]. 

The linear approximation is first reformulated as 




+ 


m 

I 


tt = 


dX w 

dp a o Pa 


<><(>„ 


(4.2.4) 


and the function <p a is chosen as 



(4.2.5) 


where r is any real number. 
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Prasad[52] also gave an alternate formulation for r = 0. The real number 
r. is a controlling parameter for the approximation. There is no obvious choice 
for r in a general problem. 


4.2.5 Generalized Hybrid Approximation 


Woo[54] combined the concepts of the conservative approximation of 
Starnes and Haftka[51] and the Generalized Inverse Power approximation and 
defined a Generalized Hybrid approximation as 



+ 



(4.2.6) 


where 


r = 


if 


g - h if 


d P a 


dX 


(*) 


G Pa 


^ 0 


< 0 


g being a real number and h being a positive integer such that 
g ^ 0 and g - h <. — 1. The choice of g and h is again not obvious, though 
the fact that larger values for g and h make the approximation more 
conservative may provide some guideline. Since this approximation too is 
only applicable to real quantities, it is applied to the real part and/or the 
imaginary part of the eigenvalue as required in particular applications. 
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4.2.6 Reduction Method(RDN) 


A different approach to derivative based approximations is based on the 
Taylor series approximation to the eigenvectors. This consists of reducing the 
original eigenvalue problem to a series of smaller order eigenproblems and is 
inspired by Noor's concept of global approximation vectors[55]. The concept 
of global approximation vectors is extended here to general matrices. 

Noor's formulation is simplified by assuming that the eigenvectors can be 
treated as linear functions of the design variables. For the sake of simplicity, 
consider only one design variable. Let 


u 


<*) 


-»?> 


+ yu (k a = [ uq ° 



(4.2.7) 


and substituting eq. (4.2.7) in eq. (4.1.1), we get, 


»»] <;> = uu?'> <>] <;> 


(4.2.8) 


If the above equation is premultiplied by [v^ vW] r , then we have 


P{y) = 


(4.2.9) 


where 


- [»f v:‘>] t a[ U w »«] 
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and 


Q - C»f> 


v (k >m*> ««] 


Eq. (4.2.9) is a 2 x 2 linear eigenvalue problem and can be solved almost 
effortlessly. Of the two eigenvalues of eq. (4.2.9), the one closest to the linear 
approximation is chosen. In the case of multiple design variables, the only 
change needed is to replace vW and uM by the respective derivatives in the 
direction of change in design. It can be proved that, if the eigenvectors are 
linear functions of the design variable, eq. (4.2.9) gives an exact eigenvalue 
when u 0 and v 0 are normalized such that their first derivatives are respectively 
orthogonal to them. This approximation will be referred to herein as the RDN 
approximation. 


4.3 Rayleigh Quotient Based Approximations 


The Rayleigh quotient for the general eigenproblem given by eqs. (4.1.1) 
and (4.1.2) is defined as 


R(x, y) 



(4.3.1) 


When all the eigenvalues of matrix A are distinct, then the Rayleigh 
quotient R(x, y) has a stationary value at x = and y = vM (the right and 
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left eigenvectors associated with the eigenvalue A.W ) for k = 1,2, ... , A/. 
Further, this stationary value is equal to XW[56]. That is, 




,1‘V 


(4.3.2) 


The approximations in this section are based on the above property. They 
seek to approximate the eigenvectors uW and vW and use them in the 
Rayleigh quotient to calculate the approximate eigenvalue. 


4.3.1 Rayleigh Quotient with Nominal Eigenvectors(RALI) 

This approximation is obtained by simply using the nominal eigenvectors 
in the Rayleigh quotient. 


K a ~ v 0 Au 0 


(4.3.3) 


Equivalently, 


4‘» = X<*> + ,f )r AAu?» 


(434) 


The second expression is cheaper to compute when AA is sparse. We will 
refer to this approximation as the RAL1 approximation. 
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4.3.2 Rayleigh Quotient with Linearly Approximated 
Eigenvectors(RAL2) 


The left and right eigenvectors are approximated from those of the nominal 
matrix using a 2-term Taylor series: 


„W _ Jk) , % ..{k) A 

U a - U 0 + 2, U a A p a 


a ~ 1 



= vW + V „W An 
Vg + L V a Ap (I 


a= 1 


(4.3.5) 


The computation of the eigenvector derivatives uM and is discussed 
in Chapters 2 and 3. 

These linearly approximated eigenvectors are then used in a Rayleigh 
quotient to generate an approximate eigenvalue. Our RAL2 approximation is 
therefore given by 



v!* )7 AuW 

A k)T ui k) 


(4.3.6) 


4.3.3 Rayleigh Quotient with Perturbed Eigenvectors(RAL3) 

In this algorithm, the perturbation AA in the matrix is used to evaluate the 
perturbations in the right and the left eigenvectors either by Adjoint or Direct 
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Method, which are analogous to the Adjoint and Direct methods described in 
Chapter 2. In the Adjoint method, assuming that the eigenvectors are 
normalized according to eq. (2.2.3) and (2.2.4), the perturbations in the right 
and left eigenvectors, 5uW and ovW respectively, are calculated as follows: 


5u' k > = £ 

7 = 1 


and 8v^ = L f k y 

7=1 


kj v o } 


where 


Vq ) 7 "aA u<*> 

(■ 4 A) - 


k *j 


vjf )7 AA uff 


k *] 


and 


e kk ~ f kk ~ 


y p ..(/) 

. ^ e kj u 0 rn 
7=1 

j*k 


(4.3.8) 


(4.3.9) 


(4.3.10) 


(4.3.11) 


McCalley[57] used this approach for error analysis of real symmetric 
eigenvalue problems. Chen and Wada[58,59] and Chen and Garba[60]used an 
equivalent formulation for real symmetric matrices and obtained eigenvalue 
approximations for all eigenvalues of the matrix simultaneously. Faddeev and 
Faddeeva[61] applied this approach to general matrices to improve the 
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accuracy of approximate eigenvalues. Meirovitch and Ryland[62] presented 
an extension to second order perturbations for general matrices. 

In the Direct method, assuming that the right eigenvectors are normalized 
according to eq. (2.2.4), the perturbations in eigenvectors are calculated as 
follows. In eq. (4.1.1), substitute eq. (4.1.3) and 

\ {k) = X ( 0 k) + 8X (k) 

u (k) = u ( 0 k) + 8u w (4.3.11) 

v W = vj*> + 5v w 

to get, after ignoring second order perturbation terms, 

A 0 5u W + AAujf* = X ( 0 k) 8u (k) + 8X (k) u ( 0 k) (4.3.12) 

Eq. (4.3.12) is identical to eq. (2.3.24) when the derivatives are replaced by 
perturbations. Hence, the same solution methods described in Section 2.3 can 
be used to solve eq. (4.3.12) for the perturbations in the eigenvector. Thus, 

5 u (k) = 0 (4.3.13) 


and 


B 0 Sy = 8r 


where 


(4.3.14) 
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B 0 - [Aq ^ 


I I 


.(*) i 

J 0 Jm-th column deleted 


5y 



with m-th element deleted 


r = - AAujj 0 (4.3.15) 

The perturbations in left eigenvectors are obtained similarly using forward 
substitution. 

The perturbations in the right and left eigenvectors are used to 
approximate the eigenvectors. Thus, 

»w - u<*> + Bu<‘> 

vf' - v«‘> + S v** 1 (4 3 16) 

The approximate eigenvectors are then used in the Rayleigh quotient as 
in eq. (4.3.6) to form an approximation to the eigenvalue. Romstad, et. al [63] 
presented both the Adjoint method approach and an equivalent Direct method 
approach of this approximation for real symmetric matrices. However, they 
obtained different numerical results with the two approaches because of an 
error in their expressions for eigenvector perturbation by Adjoint method. 
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4.3.4 Rayleigh Quotient with One-step Inverse lteration(RAL4) 


The Inverse Iteration method has been recognized as a powerful tool for 
accurate computation of eigenvectors[47,64], The unusual feature of the 
Inverse Iteration method is that accurate eigenvectors can be computed even 
when the eigenvalue is not known accurately as long as the eigenvalue is 
close enough to the correct eigenvalue. This feature can be used effectively 
to improve the accuracy of a rapid but rough approximation. In addition, a 
one-step inverse iteration is usually sufficient because most of the 
improvement in accuracy normally occurs in the first step and in a 
modification problem, the nominal eigenvector is available and provides an 
excellent initial iterate. 

In this algorithm, a first approximation to the eigenvalue is formed 
as a Rayleigh Quotient (RAL1) given by eq. (4.3.3). This approximation is then 
used in a one-step inverse iteration scheme to obtain approximate 
eigenvectors as follows: 


u«*> = (A - 4‘,'ir 1 
v<‘> - <A r - Or 1 vf?' 


(4.3.17) 


The approximate eigenvectors are then used in the Rayleigh quotient as 
in eq. (4.3.6) to form an approximation to the eigenvalue. We will refer to this 
approximation as RAL4. 
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Note that the evaluation of the and v(f) requires only one matrix 
factorization, since the second part of eq. (4.3.17) can be solved by forward 
substitution using the same factored matrix used in solving the first part. 


4.4 Trace-Theorem Based Approximations 


These approximations are based on well known iterative methods for 
finding the roots of a polynomial. We apply these methods to the 
characteristic polynomial, p(X) of matrix A using only one step of the iteration 
for the approximation. 

The remarkable feature of the approximations in this section lies in the fact 
that the coefficients of the characteristic polynomial need not be calculated 
explicitly. This is achieved by employing the following result, known as the 
Trace Theorem[56], 

Trace Theorem: If p(X) # 0, then 


m 


p'(V) 

pm 


= — [ Trace of (A -XI) ^ 


(4.4.1) 


where p'(X) = 


dp 

dX 
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4.4.2 Refined One-step Newton-Raphson lteration(NRT2) 

In this approximation, we utilize the second derivative of p(X) to obtain a 
better approximation. To refine the Newton-Raphson Iteration using the 
second derivative, consider the truncated Taylor series expansion 

0 = P(4 k ) = Pfral) + P'^ ( ah^a ] ~ + OSp"^)^ " ^1 ) ) 2 (4.4.4) 
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Use eq. (4.4.2) to eliminate X^ in the third term on the right hand side, to get 


0 - + MW - 

2ip'(4i)r 


(4.4.5) 


so that the refined Newton-Raphson formula is 


*(*) _>(*)_ 
A a ~ A a 1 


P(^al) P'UL*V(*i*1> 


p'( 4 i) 


2[pv3f 


(4.4.6) 


From the definition of f(X), we have 


p"(A) - P(>v)lr(A) + f 2 (A)j 


(4.4.7) 


Now, differentiating eq. (4.4.1) with respect to X, while noting that 


(A — X\) 1 = C(A — XI) 1 ] 2 
c/A 


(4.4.8) 


we get, 


f(A) = - Trace of [(A - Al) ^ 


(4.4.9) 


This expression is then used in the refined Newton-Raphson formula of eq. 
(4.4.6) to obtain our next approximation given below. 


.(» .w n4* ?> + 3/Vgf) 

a al o / ^ \ 

2f 3 (X^) 


(4.4.10) 
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Note that, in evaluating f(X) using eq. (4.4.9), the matrix multiplication need 
not be performed completely, since only the diagonal elements of the matrix 
product are needed. 

This approximation will be referred to herein as the NRT2 approximation. 

4.4.3 One-step Laguerre Iteration(LIT) 


Laguerre iteration is often used to compute eigenvalues and is known to 
have excellent convergence properties[47]. For our purposes, we need only 
one step of the Laguerre iteration. The one-step Laguerre iteration consists 
of 


^ a ^al 


r>P(K l) 


pUai) ± [(n - 1) 2 p' 2 (X a1 ) - n(n - 1)p(X a1 )p"(A a1 )] 


1/2 


(4.4.11) 


This approximation also utilizes the second derivative of p(X) . The 
derivation of eq. (4.4.11) is given in Wilkinson[47]. To obtain our 
approximation, we rewrite the above, using the trace theorem, as 


X a = X a1 - (4.4.12) 

f(X a i) ± [ - (n ~ 1 )f 2 (k a1 ) - n(n - 1)f(X a1 )] 1/2 

The evaluation of f(X) is described in Section 4.4.2. The sign in the 
denominator is chosen so as to make the denominator have the greater 
absolute value. We will refer to this approximation as the LIT approximation. 
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4.5 Other Approximations 


An alternative approach to approximate eigenvalues is taken by Paipetis 
and Croustalis[65] who developed an algorithm to approximate the coefficients 
of the characteristic polynomial which is then solved to obtain approximate 
eigenvalues. In addition to the numerical difficulties associated with the 
evaluation of the coefficients of the characteristic polynomial and its solution, 
this method severely restricts the characteristics of the system matrix. 

We will consider one more approximation called [1,1]Pade approximation. 

4.5.1 [1,1] Pade Approximation(PADI) 

We derive the [1,1] Pade approximation by geometrical construction along 
the lines of Johnson[66], Let us for the moment assume that the eigenwiui; 
to be approximated is real. 

Let ar, d be the first and second approximations to the eigenvalue 
XW. The information contained in these approximations is exploited by 
Aitken's method[67] to obtain a hopefully better approximation. We form the 
differences 


sn (k) _ (k) _ (k) 

OA a i — A a1 /V 0 

Sft(k) = >(k) _ i (k) 

0A a 2 ^a2 1 


( 4 . 5 . 1 ) 
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If XtfK X^ and X^ is a converging series, we will have 

SX a2 < SX a1 (4.5.2) 

We now draw a straight line through the points (Xf^), 8X^) and (Xi$, 5X^) 
to extrapolate to 8XW = 0 on a (XW, SXM) plot. This is illustrated in Figure 3. 
It is expected that the extrapolated value X^) would be a better approximation 
than either X^ or X^. 

The result is 



+ 


s4‘i’ - 


(4.5.3) 


Although the motivation applied only to real eigenvalues, this result can 
be immediately extended to complex eigenvalues. Using eqs.(4.5.1), eq. (4.5.3) 
is rewritten as 



nW)2 _ 

lA. a l ) Aq A a2 

^A. a i Aq A a2 


(4.5.4) 


If the approximations X^,X(^ are the linear and the quadratic 
approximations, the approximation of eq. (4.5.4) can be recognized to be in the 
form of the [1,1] Pade approximant. In the following, this is assumed and this 
approximation will be referred to as PAD1 approximation. 
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Chapter 5 

Accuracy and Efficiency of Eigenvalue 
Approximations 


5.1 Introduction 


In Chapter 4, we listed several approximations for the eigenvalues of 
general matrices. For a given application, the selection of the appropriate 
approximation usually depends on the saving of computational time that a 
given approximation entails. In the task of selecting a good approximation, 
information about the accuracy and the efficiency of computation of the 
approximations is essential. Accuracy and efficiency are not independent 
elements in the selection of an approximation algorithm. Poor accuracy 
usually translates into low efficiency in the global process. 
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In this Chapter, the accuracy and efficiency considerations relating to the 
approximations listed in the last chapter are discussed. The accuracy 
considerations are treated in Section 5.2 and the efficiency considerations in 
Section 5.3. The Conservative, Generalized Inverse Power and the 
Generalized Hybrid approximations are not studied as their accuracy is 
problem dependent and a general assessment is not feasible. 

Some of the approximations discussed have been applied to symmetric 
matrices by researchers in structural dynamics. However, there exists no 
systematic comparison of accuracy and efficiency. To the best of the author's 
knowledge, the trace theorem based algorithms have never been used for 
approximating eigenvalues of modified systems in the structural dynamics 
literature. 


5.2 Accuracy Considerations in Approximating 
Eigenvalues 


For simplicity of notation, we consider a single design variable. It will be 
obvious, however, that the results of the error analysis are applicable to the 
case of multiple design variables as well. 
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5.2.1 Order of an Approximation 


For a useful comparison of the various approximations, we define the 
order of an approximation as follows: 

Definition: If an approximation X a to an eigenvalue X is such that the error 

(X a - X) = 0[(Ap a ) s+1 ] or (k a - X)= C(Ap a ) s+ 1 for A p u - 0 (5.2.1) 

then that approximation is said to be of s-th order. 

Note that, for a rigorous estimate of the error in an approximation, 
information about both the order of the approximation as well as the 
proportionality constant C in eq. (5.2.1) is needed. However, the order of the 
approximation is usually the most important property of the approximation and 
the proportionality constant is useful only when comparing approximations of 
the same order. In this work, attention is focused on the order of the 
approximation and the proportionality constant is discussed only in examples. 

5.2.2 First Order Approximations 

Among the approximations described in the last chapter, the linear 
approximation(LIN) and the Rayleigh quotient with nominal eigenvectors 
(RAL1) are first order approximations. 
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It is easy to show that the Linear approximation is of first order. Consider 


the two term Taylor series for A.W with remainder given by 


p a + 4 ^- hww pf 


AW 


°Pa 


d 2 p. 


(5.2.2) 


where P a £ C a ^ P a + A Pa • 

Comparing this to the linear approximation, eq. (4.2.1), we have the error 

4 k) - = °( A P«) (5.2.3) 


Hence the linear approximation is a first order approximation. 

To find the error in the RAL1 approximation, subtract eq. (4.3.4) from the 
exact expression of eq. (4.1.10). Ignoring the third order terms, we have 

»[f» r AAAu« k > + Av Wr A 0 Au« + A»W r AAu<f> 
vjf^AAu^O^W* 1 + Av^uj,* 1 ) - X 0 Av W W*> 

1 + Av^' 7 Uq + VgAu + Av^^Au^ 

Considering that, to the first order, 

AA, Au, Av = 0(Ap a ) 

we have 

X^ — = 0(Ap a ) (5.2.5) 


X 


W 


X(k\ 


PALI 
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Thus, RAL1 is a first order approximation. 


5.2.3 Second Order Approximations 

The Quadratic approximation and its improvement by [1,1] Pade 
approximation(PADI) are the second order approximations we described. To 
show that the quadratic approximation is of second order, consider the three 
term Taylor series with remainder given by 


*<*> - 4' + P„ + 

u rci 


1 ... .2 . 1 


3,W 


2 , 2 
°Pn 


(Ap (1 ) + r (C a )(Ap u )° (5.2.6) 


Comparing this to eq. (4.2.2), we have the error in the quadratic 
approximation as 


K a k Q\JAD u ^PJ 


(5.2.7) 


showing that the quadratic approximation is a second order approximation. 
The PAD1 approximation is an improvement on the quadratic approximation 
and so it is at least of second order. 

The linear, quadratic and the PAD1 approximations are applicable to all 
functions and do not take advantage of the special properties of the eigenvalue 
problem. All the other approximations we are going to discuss are developed 
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specifically for approximating the eigenvalues and it will be shown that they 
achieve higher accuracy with less computational effort. 


5.2.4 Third Order Approximations 


The reduction method RDN, the Rayleigh quotient based methods RAL2 
and RAL3 and the one-step Newton-Raphson iteration NRT1 are the third 
order methods of approximation that we considered. As the first three 
approximations RDN, RAL2 and RAL3 are closely related, we will derive the 
order of only the RAL2 approximation and infer the orders of the RAL3 and the 
RDN from this derivation. 

Recall that in the RAL2 algorithm, the eigenvalue is approximated by using 
linearly approximated left and right eigenvectors in the Rayleigh quotient. 
Hence, we may write the two term Taylor series with remainder for the 
eigenvectors as 


= „?> 
•<*> = v 


+ Yv!l(C„)(Ap a ) 2 


(5.2.8) 


Hence, the approximate eigenvectors can be rewritten in terms of the exact 
eigenvectors as 
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(5.2.9) 


„f - » w - 4y*i<y<W 


,W = „(*) - JL„M 


V 

2 , «« 


(y(4P„r 


Substituting these expressions in eq. (4.3.6) and using the eqs. (4. 1.1-2), we 
get 



vj, k,T Au ( a k> 

-fer 


X W (1 - Z W ) + 0(Ap*) 
1 - Z^ + 0(Ap<) 


(5.2.10) 


where 


-W 


- + ,<‘>V*> a <y]<Ap„, : 


(5.2.11) 


Carrying out the long division and putting Xif) = X$b L2 , we find the error 
in the RAL2 approximation as 


jW jW — DfAn^l 

A ^RAL2 U ( A P«) 


(5.2.12) 


It is hence proved that the RAL2 approximation is of third order. 

It may be recalled that in the RAL3 approximation, we approximate the 
eigenvalue by a Rayleigh quotient using left and right eigenvectors that were 
obtained by using first order perturbations whereas in the RAL2 algorithm, we 
approximated the eigenvectors using their first derivatives in a two-term 
Taylor series. There is mathematically no difference between these two 
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methods and their difference lies only in the computational algorithms. It 
follows then that RAL2 is also a third order approximation. There is indeed 
little difference in the accuracy of the two approximations when they were 
tested on example matrices. 

The order of the reduction method approximation(RDN) is difficult to 
obtain algebraically. However, we can infer the order of the RDN 
approximation by comparing it to the RAL2 approximation. We first note that 
the chief characteristic of the RAL2 approximation is that the eigenvector of 
the modified matrix is assumed to be in the subspace consisting only of the 
linearly approximated eigenvector. The RDN approximation is more flexible 
in that the eigenvector of the modified matrix is assumed to be in the subspace 
consisting of the original eigenvector and its derivative. Thus, the RDN 
approximation can be expected to be somewhat better than the RAL2 
approximation. Hence, the RDN and the RAL2 approximations are expected 
to be of the same order but possess a different proportionality constant in the 
sense of eq. (5.2.1). This conclusion is validated by several numerical 
experiments. 

To derive the order of the one-step Newton-Raphson iteration(NRTI), we 
first write the Taylor's series for the characteristic polynomial, p(X) and its 
derivative, p'(k), as 

P^ai) = P(^ W ) + (4*? ~ ^ ( V + 0.5(4*/ - X W ) 2 p" + - (5.2.13) 

and 
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p'(4*?) = p' + <4*i - *■ {k) )p " + o- 5 (4 k i - ^ k) ) 2 p"’ + - 


(5.2.14) 


where primes denote derivatives evaluated at XW. Then the Newton-Raphson 
Algorithm of eq. (4.4.2), modified by by subtracting XM from both sides and 
noting that p(X M) = 0, 


»(*) _ i(*) = _ *(*) 

/la A A q 4 A 


al 

(XjJ? - X W )p' + 0.5(Xj$ - X W ) 2 p" + 
p' + (X$ - X W )p" + - 


(5.2.15) 


may be written as 


x< fc) - x<*> = 


0.5(X^ 1 ) - X W ) 2 p" + (1/3)(XJJ - X W )V" + 
p’ + (X<*> - X W )P" + ... 


(5.2.16) 


giving 

XW-XW-(XW-XW) 2 ^- (5.2.17) 

We have chosen the initial approximation X^ to be the RAL1 approximation, 
which has already been shown to be a first order approximation. Hence, 

(4*1 - *•' W ) = 0(Ap 2 ) (5.2.18) 

so that 
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OjRn - *“>) 


(5.2.19) 


- O(APa) 

establishing that the NRT1 approximation is a third order approximation. 


5.2.5 Higher Order Approximations 


The remaining approximations, RAL4, NRT2 and LIT are all fifth order 
approximations. We proceed to obtain the order of these approximations. 

To get the order of the RAL4 approximation, we follow Ostrowski's 
approach[68] using, however, the simplifying assumption that all the 
eigenvalues are well-separated. Let U and V denote the matrices whose 
columns are the right eigenvectors u and the left eigenvectors v respectively 
of the matrix A. Let A denote the diagonal matrix of eigenvalues. From the 
biorthogonal property of the left and right eigenvectors normalized as given 
by eq. (2.2.3), we have 

V r AU = A and V 7 U = I (5.2.20) 


Define 


n (k) _ V T Ak) 
% ~ v u 0 • 






i" 1 - v’uf. 


F (k) 


U r vf 


(5.2.21) 


From eqs. (5.2.20), we have 
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V r (A - Xjfll)U = A - 


(5.2.22) 


so that 


(A = U(A - 


(5.2.23) 


Using the above and the definitions of eq. (5.2.21) after premultiplying the 
first part of eq. (4.3.17) by V r , we obtain the relation 


tip) - (A - Xjflir 1 up) 


(5.2.24) 


From the second part of eq. (4.3.17), we can obtain, in a similar manner, 


tf) - (a - 4VT 1 tf) 


(5.2.25) 


Using eqs. (5.2.20-21), the RAL4 approximation can now be written as 


Uk) , 


(5.2.26) 


Substituting the expressions (5.2.24) and (5.2.25) in eq. (5.2.26), we have 


£ 

/=1 (A. (;) - 4V) 2 

£ H 

/=1 (a. (/) - 4V) 2 


(5.2.27) 
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Subtracting the exact eigenvalue A.W from both sides of this expression, we 
get after considerable algebra, 


- x'"' 

(X'*» - X <‘>) 2 


n (k® 
v 

/'= 1 

i*k 




(X« - o 


<*K2 


«M? + .2 nW 


i- 1 
i*k 


(X (t > - xj , 1 ?) 2 
<A W - 4\) 2 


(5.2.28) 


Then, if the eigenvalues are weii separated and if }J$ is a reasonably close 
approximation of the eigenvalue , the second term in the denominator will 
be negligible compared to the first. So, we may write 


k {k) - \ {k) 

oi - * w ) 2 


n 

I 

/'= 1 
/ ^ A 



e (*)«(*) 
bOA ^OA 


(5.2.29) 


However, is taken to be the RAL1 approximation. Hence we have, 




(5.2.30) 


Subtracting from both sides as before, 
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4*, 1 - x '*» = 


£ - xW 

/ = 1 

i*k 

eM? + SW 

/= 1 

i*k 


(5.2.31) 


Or, 




£ (x« - x«kM> 

/=i 

/>* 

(Xg? - x<*>) 


- £ *M> 
/ = 1 

i*k 


(5.2.32) 


When the eigenvalues are well-separated and X^ is close to the exact 
eigenvalue A.W, the first term is dominant, so that 


n 

I (X 
/=1 

/*/c 


(/) 




.<*)_<*) ~ 
SO/c^O/c = 




w - X W) 


(5.2.33) 


Substituting eq. (5.2.33) in eq. (5.2.29), 


80 


xf> - X<*> 

< x £? - x <*>) 2 


g (x w - 

(X» - xj$) J 

/** 


cxffl 


- X<*>)- 


£ (X» 

/•=1 

i*k 


X^KiM 


Now, putting ^ = ^$Ai.4 and = ^$Ali> we have 


A rt>4L4 



(5.2.34) 


(5.2.35) 


Since RAL1 is a first order approximation as shown in Section 5.2.2, 

PL-l" 1 ! * O(^Pa) I 5236 ’ 


proving that the RAL4 approximation is of fifth order. 

The derivation of the order of the NRT2 approximation is analogous to that 
of NRT1 approximation given in the last section if eq. (4.4.10) is used in place 
of eq. (4.4.2). After considerable algebra, we get 

Xj&„ - X» = 0[(X$ - X (k) ) 3 ] (5.2.37) 

For the Laguerre Method, Wi!kinson[47] gives 

X?,V - X ( *> - 0[(XS> - X w f] (5.2.38) 

Since the initial approximation used in both these algorithms is of first order, 
we have 
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( 5 . 2 . 39 ) 


&NRT2 ~ ~ O(APq) 

(X[5V - X (k) ) = O(ApS) (5.2.40) 

so that both the NRT2 and the LIT approximations are fifth order 
approximations as stated. 


5.2.6 Validation of the Theoretical Results 

Taking logarithms of both sides of eq. (5.2.1), we have 
log (X a - X) - (s + 1) log (Ap a ) + constant (5.2.41) 

so that if the error in an s-th order approximation is plotted against the change 
in design variable on log-log scale, one must obtain a straight line with slope 
(s + 1). 

The orders of approximations obtained in the last section are verified by 
numerical experiments, summarized in the following figures. In order to 
minimize the effect of round-off errors as much as possible, a small matrix of 
order 5 is used for validation of theoretical results. The matrix elements are 
the quadratic polynomials of a design variable where the coefficients are 
generated using a random number generator. The errors in the approximate 
eigenvalues are in comparison to the exact eigenvalue that is obtained by the 
QR algorithm and improved by using the eigenvectors computed by inverse 
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Figure 5. Relative Percentage errors in the Real Part of the Eigenvalue 5x5 Matrix 
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iteration in a Rayleigh quotient. Figures 4, 5 and 6 show the errors in the 
absolute value, the real part and in the imaginary part of an eigenvalue of the 
5x5 matrix plotted on a log-log scale against the change in design variable. 
In all the following figures, 10' is represented in the FORTRAN notation IE/. 
The slopes of the straight line segments in Figure 4, 5 and 6 agree very closely 
with the orders of the approximations derived in Section 5.2. The difference 
in accuracy between the high order and the low order approximations is 
clearly apparent. 

Approximations are also applied to one of the eigenvalues of a larger 
matrix of order 40, generated in the same manner as the smaller matrix of 
order 5 before. The results obtained are shown in Figures 7, 8 and 9. The 
results for a flutter analysis matrix are shown in Figures 10, 11 and 12. The 
generation of the flutter analysis matrix is described in Section 3.2. The design 
variable used in the Figures 10, 11 and 12 is the reduced frequency, 
k defined as where a> is the vibration frequency, b the blade 

semi-chord and V the air speed in far field. 

The deviations from straight line behavior appearing in Figures 10, 11 and 
12 at small values of A p a are typical instances of round-off error with which 
we are not concerned. 

In all these Figures, approximations of equal order are indicated by 
straight lines of equal slope. The proportionality constant, which reflects 
accuracy, is indicated by the vertical position of the corresponding straight 
line. 
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RELATIVE ERROR IN IMAGINARY PART 
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It is seen that both the first order approximations have about the same 
accuracy and that the improvement achieved by the PAD1 approximation over 
the quadratic approximation is marginal. The reduction method approximation 
with good consistency shows substantially more accuracy than the other third 
order methods. Among the fifth order approximations, the NRT1 
approximation is consistently poorer in accuracy than others. 


5.3 Efficiency Considerations in Approximating 
Eigenvalues 


Efficiency is probably the most important consideration in the comparative 
evaluation of the various approximations, particularly in the context of design 
optimization. The comparison is once again made in terms of variables which 
significantly influence the cost of computing an approximation. In addition to 
the size of the matrix n, the number of design parameters m and the number 
of eigenvalues of interest / that we considered in Chapter 3, we have an 
additional variable in the approximation context and this is the number of 
design points d at which an approximation is sought based on the same 
nominal design. 

Tables 4 and 5 present the operation counts for all the approximations 
studied. The operation counts include the necessary computations of the left 
and right eigenvectors at the nominal design represented by the matrix A 0 , if 
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Table 4. Operation Counts for First and Second Order Approximations 


Approximation 

LIN 

RAL1 

Approximation 

QUAD 


First Order Approximations 

Operation Count 
l{-^-n 2 + k mn 2 + dn) 
l{-^-n 2 + k dn 2 ) 

Second Order Approximations 

Operation Count 

/—- + / ) n 2 k + lmn 2 { 2k + 1) + din 2 

- Direct-Adjoint Method 

-j-J n 3 + (k + 1)mn 3 + / n 2 k + din 2 

- Adjoint Method 


PAD1 


Same as QUAD 


Table 5. Operation Counts for Third and Higher Order Approximations 



Third Order Approximations 


Approximation 

Operation Count 


RDN 

/[-^- — f 2 (k + 1 )mn 2 + 2dn 2 ] 

O 

Direct Method 


-^-n 3 + Ip 2 [{k + 2 )m + 2d] 

Adjoint Method 

RAL2 

/[-^- — t- 2 (k + 1)mn 2 + dn 2 ] 

Direct Method 


-^■n 2 + ln 2 [{ k + 2)m + d] 

Adjoint Method 

RAL3 

d/[-^- + 2 (k + 1)n2] 

Direct Method 


-^■n 2 + dln 2 (K + 3) 

Adjoint Method 

NRT1 

din 2 



Higher Order Approximations 


Approximation 

Operation Count 


RAL4 

n 2 


NRT2 

din 2 


LIT 

din 2 
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they are significant. But they do not include the operations for the calculation 
of the nominal eigenvalues since these do not affect the comparison of the 
efficiency of the different approximations. When the derivatives of 
eigenvectors or the second derivatives of the eigenvalues are needed for an 
approximation, operations counts for both the Adjoint method algorithm and 
Direct method algorithm are given separately. The details of the Adjoint and 
the Direct methods are discussed in Chapters 2 and 3. 

Note that the computational expense of RAL1, RAL3, RAL4 and all the fifth 
order approximations is independent of the number of design variables. 


5.4 Discussion of Approximations 

5.4.1 Case When No Derivatives are Available 

The results depicted in Figures 4-12 show that both the first order 
approximations we studied, LIN and RAL1, are practically identical in 
accuracy. Comparing their operation counts, we conclude that, when first 
order accuracy is acceptable, LIN approximation is preferable if the number 
of design variables is small and the number of design points for approximation 
is large. When the number of design variables is large, the RAL1 
approximation is more efficient. 
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The experience of numerical experiments also shows that the improvement 
achieved by the PAD1 approximation over the quadratic approximation is 
marginal so that we have two second order approximations of nearly the same 
accuracy. The behavior of the [1,1] Pade approximation based on initial 
approximations other than the linear and the quadratic has not been studied. 
However, the operation counts show that the evaluation of the second 
derivatives is very expensive. The third and higher order approximations are 
not only more accurate but are also more efficient so that the second order 
approximations QUAD and PAD1 can be completely dropped from 
consideration. Among the third order methods, RAL2, RAL3 and NRT1 have 
similar or same accuracy while the Reduction method RDN shows higher 
accuracy in most cases, sometimes close to that of some of the fifth order 
approximations. This is explained by the fact that the reduction method 
approximates the eigenvectors in a subspace of two vectors whereas RAL2 
and RAL3 approximations use a subspace of only one vector. However, the 
reduction method is also more expensive than RAL2 and RAL3 as shown by 
the operation counts in Table 5. 

Hence, among the third order methods, the trade-off between accuracy 
and efficiency determines the choice of the approximation. When accuracy is 
more important than efficiency, the reduction method is chosen over the 
others. When efficiency is more important, we have the choice between RAL2 
and RAL3. The Newton-Raphson approximation(NRTI) is dropped from 
consideration because there are higher order methods which give higher 
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accuracy with the same computational expense. Note that the Direct method 
computation of the RAL3 approximation is more expensive than the Adjoint 
method computation except for very few eigenvalues of interest and very few 
design points of approximation. The choice between RAL2 and RAL3 is similar 
to that between LIN and RAL1 in the first order case. RAL2 approximation is 
preferable if the number of design variables is small and the number of design 
points for approximation is large. When the number of design variables is 
large, the RAL3 approximation is more efficient. 

The higher order approximations are particularly efficient when the 
number of design variables is large and the number of design points and the 
number of eigenvalues of interest is small. The RAL4 algorithm is somewhat 
more efficient than the other two. Among these approximations, for the cases 
tested, the RAL4 and the LIT methods give better accuracy than NRT2. As the 
NRT2 approximation is no cheaper than the LIT approximation, this eliminates 
the NRT2 approximation from consideration. Whether there is any 
considerable difference in accuracy between the RAL4 and the LIT methods 
requires further investigation. In the absence of any such difference, the RAL4 
approximation may be considered to be the best higher order approximation 
available in terms of accuracy and efficiency. 

The computational cost of the higher order approximations escalates 
rapidly to equal that of the exact computation of eigenvalues. Comparing the 
operation count of the RAL4 approximation and the exact computation, we 
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note that the RAL4 approximation is more expensive than the exact 
computation when the product dl is above 30. 

5.4.2 Case When Derivatives are Available 

Most design optimization algorithms require first derivatives of 
constraints. Design optimization of dynamic response in structures almost 
always involves constraints on the eigenvalues and sometimes constraints on 
eigenvectors are also involved. In such cases, the first derivatives of the 
eigenvalues and perhaps eigenvectors are already available free for use in 
approximations so that the operation counts for derivative-based 
approximations given in Tables 4 and 5 are reduced, affecting some of the 
conclusions. In this section, we discuss the relative merits of the 
approximations when the first derivatives are already available. 

We first consider the case when constraints are placed only on the 
eigenvalues so that the first derivatives of only eigenvalues are available free. 
In such a case, the operation count for only the linear approximation is 
affected and is shown in Table 6. The linear approximation is now an 0(n) 
process and is practically free in terms of computational expense compared to 
any other approximation. Thus, when constraints are placed on the 
derivatives of eigenvalues, the linear approximation is the best approximation 
unless particularly high accuracy is required. If high accuracy is desired, the 
conclusions of Section 5.3.1 still hold. 
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Table 6. Operation Counts for First Order Approximations when Eigenvalue Derivatives are Free 


Approximation 

LIN 

RAL1 


First Order Approximations 

Operation Count 
Idi 7 

l(-^-n 2 + Kdn 2 ) 
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We next consider the case when the constraints are placed on both the 
eigenvalues and eigenvectors so that all first derivatives are available free. In 
such a case, the operation counts for all derivative-based approximations are 
affected and are shown in Tables 7 and 8. 

The derivative based approximations except the quadratic approximation 
are much more attractive when all the first derivatives are available. The 
linear approximation is again practically free so that it always makes sense to 
use linear approximation in terms of efficiency. The third order 
approximations RDN and RAL2 are now 0(n 2 ) processes and are substantially 
cheaper than the fifth order approximations which are still 0(n 3 ) processes. 
However, the RDN approximation is now twice as expensive as the RAL2 
approximation and the additional computational expense of the RDN 
approximation is less easily justified even though it is somewhat more 
accurate. The quadratic approximation is again more expensive than the more 
accurate third order approximations. Operation counts for the fifth order 
approximations are not affected by the availability of the first derivatives and 
hence the conclusions regarding the same are also unaffected. 
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Table 7. Operation Counts for First and Second Order Approximations when All First Derivatives 
are Free 


First Order Approximations 


Approximation 

LIN 

RAL1 


Approximation 

QUAD 


PAD1 


Operation Count 

Idn 

Idi^-n 2 + Kn 2 ) 

Second Order Approximations 

Operation Count 

I n 2 k + lmn 2 K + din 2 

- Direct-Adjoint Method 

(k + 1 )mn 2 + / ( 2 *) n 2 K + din 2 

- Adjoint Method 
Same as QUAD 
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Table 8. Operation Counts for Third and Higher Order Approximations when All First Derivatives 
are Free 


Approximation 

Third Order Approximations 
Operation Count 


RDN 

4dln 2 


RAL2 

din 2 


RAL3 

~n 3 -f d/n 2 ( k + 3) 

Adjoint Method 

NRT1 

din 2 


Approximation 

Higher Order Approximations 

Operation Count 


RAL4 

n 3 

dl <ir> 


NRT2 

din 3 


LIT 

din 3 
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Chapter 6 
Conclusions 


The large computational expense associated with the flutter optimization 
of a cascade of rotating blades motivated this study. The problem of 
computational expense is attacked from two fronts, sensitivity analysis and 
approximations. The existing literature was surveyed in both fields and 
improvements are suggested. General recommendations for the selection of 
the most efficient algorithms are presented. 

The normalization of the eigenvector needs to be properly related to its 
derivative. In practice, this means that the derivative of the eigenvector is to 
be normalized before it is used, to conform to the normalization of the 
eigenvector itself. When the eigenvector is not normalized in a unique 
manner, its derivative cannot be evaluated. It has been shown that fixing one 
of the components of the eigenvector is the best normalizing condition for 
computation of the derivative. 
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In the sensitivity analysis part, the algorithms presently available for 
computing exactly the derivatives of eigenvalues and eigenvectors are 
classified into Adjoint and Direct Methods. Adjoint Methods use both the left 
and the right eigenvectors whereas the Direct Methods use only the right 
eigenvectors. The Adjoint Methods and the Direct Methods found in the 
literature are extended to apply to eigenvectors normalized in the manner 
described above. Algorithms that compute approximate derivatives are not 
studied as their implementation is problem dependent or complicated. 

The Adjoint and the Direct methods are examined for their efficiency under 
different sets of conditions. The choice reflects whether the solution of the 
adjoint problem is worth the extra computational expense. The solution of the 
adjoint problem is shown to be cheaper than the solution of the direct problem 
because left and right eigenvectors can be calculated using the same factored 
matrix. It is concluded that if only the first derivatives of eigenvalues are 
required, the solution of the adjoint problem is worth the expense since the 
Adjoint method is superior to the Direct Method. When first derivatives of 
eigenvectors are also required, the decision is dependent on the problem size, 
the number of design variables and the number of eigenvalues of interest. The 
Direct method is more competitive if the number of design variables is large 
and the eigenvalues of interest are few. When the first and second derivatives 
of eigenvalues are required, similar considerations hold. It is also shown that 
once the first derivatives of eigenvectors are calculated, the second 
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derivatives of eigenvalues are calculated more efficiently by the Adjoint 
method than by the Direct method. 

In the approximation part, many existing approximation methods are 
applied to general matrices. Some new approximation methods, inspired by 
the computational techniques in linear algebra, are proposed. Noor's concept 
of global approximation vectors is simplified and then extended to general 
matrices to arrive at another approximation called the reduction method. 

The approximations are classified as Derivative Based, Rayleigh Quotient 
Based, Trace Theorem Based and [1,1]Pade approximations according to their 
theoretical origin. In terms of accuracy, the approximations are reclassified 
according to the order of the error expected from the approximation. The 
approximation methods are also examined for computational expense as this 
is a significant issue in the selection of a particular method. At each order of 
accuracy, the approximations are compared in terms of their efficiency and 
general recommendations are made. Additional recommendations are made 
for the case when the derivatives are already available. In particular, it is 
concluded that the quadratic approximation is inferior to many other 
approximations both in accuracy and efficiency. 

The analysis performed in this work is applicable also to the generalized 
eigenvalue problem in a straightforward manner so that all the conclusions are 
valid for the usual structural stability problem. However, the conclusions are 
limited by the assumption of distinct and well-separated eigenvalues of 
interest. The sensitivity analysis and approximation methods for multiple or 
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closely spaced eigenvalues is fraught with difficulties, numerical as well as 
theoretical. The multiple eigenvalue case is suggested as a topic for further 
research. 
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